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DERIVE & CAS-TI User Group Meeting 2004 

User Group Meeting 2004 was announced for Saturday, 17 July, during TIME 2004 Conference in 
Montreal, Quebec, Canada. 

DUG-members who attended the conference followed the announcement and joined the meeting. We 
spent two interesting and inspiring hours together. 

The editor and president of the DUG gave a report about the activities of the DUG so far: 

The assembly accepted and appreciated the editor's proposal to have no membership dues for 2004 
and the following years due to the fact, that the newsletter will not be printed and shipped by snail 
mail, but can be downloaded from the DUG-homepage. 

Electronic form of the DUG publications meets the wishes and suggestions of many members. About 
50 new members joined the User Group in 2004. Welcome to them all. 


Since founding the DUG in 1991 we have published 54 newsletters 
containing 

more than 2200 pages 
386 contributions submitted by 
149 authors from 
27 countries. 

We had 592 requests and answers in the User Forums and presented 

239 CAS-related books on the book shelf. 

Among others Johann Wiesenbauer provided 27 Titbits from Algebra and Number Theory 
which could fill a very special book on this topic. We were very happy that we could express our 
gratefulness to Johann personally at the meeting. 

We found that the Derive-Users have been very productive in delivering papers, the CAS-TI-Users are 
hesitating - with a few exceptions. So we have to encourage this group to go to public with their find- 
ings in the future. 

Ideas and intentions for the future are: 

♦♦♦ Continuing publication of revised versions of Newsletters from 1991 - This service 

is very much appreciated by the members and is an exciting task for the editor. 

❖ Providing an extended index containing all articles and contributors. Later we would like 
to add short abstracts (2 sentences) and links between related contributions. 


(continued on page 3) 
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Dear DUG-Members, 

First of all I'd wish to welcome our 50 new mem- 
bers who joined DUG during 2004. We hope that 
you will find many information and inspiring ideas 
in our publications and we would be very pleased if 
some of you would submit one or the other contri- 
bution for publication. 

As I reported at the DUG-Meeting in summer we 
would like to encourage our handheld users to 
share their experiences with us. I had many positive 
reactions on the fact that the Newsletter can be 
downloaded from our website and that I revise the 
old newsletters from our first years. I can say that 
this is a very interesting task to recognize how 
things have 
changed within a 
few years. I am 
sure that we will 
discover many 
new qualities of 
the many contri- 
butions from the 
early nineties. I 
was relly moved 
by revising 
DNL#4, when I 
found again a Call for Papers for our first DERIVE 
Spring School 1992 in Krems. Those were the days 
when it all has begun. In 1991 the Austrian gov- 
ernment purchased the nationwide DERIVE licence 
for our secondary schools and we are very happy 
that some weeks ago our government renewed the 
licence contract for DERIVE 6.10. 


Bjoem's article was written before times of De- 
rive's slider bar and I tried to include this valueable 
tool. I left his article in Danish and added some 
English comments. 

Just recently I 
received a mail 
from our very 
productive 
member Don 
Phillips who 
sent a very 
impressive 
paper on "Two 
stage least 
squares", an- 
other mail from Josef Lechner containing two arti- 
cels for publication and one mail from Canada 
announcing a paper about "Diophant-ive Polyno- 
mials". Lorenz Kopp did not only send his fine tool 
for generating tree diagrams but also a bundle of 
simulations for random experiments with very 
interesting graphic representations. 

Please take also notice of the extended report on the 
DUG-meeting 2004 in Montreal. 

I'll take the occasion to thank 
Walter Wegscheider and Ben- 
jamin Kaineder for putting the 
DNLs and other files on the 
website. 

Finally I wish you and your families a Merry 
Christmas Time and a Happy, Healthy and Success- 
ful New Year 2005. 





This issue contains two fine articles from Giora & 
Nurit (Quadratic Approximation) and Karsten 
Schmidt who fulfills his promise to present exam- 
ples for application of Moore-Penrose Matrices. 
We have a Danish contribution which deals with 
Quadratic Optimization and shows that this topic 
can be treated even in secondary school. 


With best regards from my grandchildren and from 
Noor and Josef 



Download all DNL-DERIVE- and Tl-flles from 

http : //www . austromath . ac . at /dug/ 
http : / /www . bk- teachware . com/ main . asp?session=37505 9 
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The DERIVE-NEWSLETTER is the Bulle- 
tin of the DERIVE & CAS-77 User Group. 
It is published at least four times a year 
with a contents of 44 pages minimum. The 
goals of the DNL are to enable the ex- 
change of experiences made with DERIVE 
and the 77-89 /92/ Titanium! Voyage 200 as 
well as to create a group to discuss the 
possibilities of new methodical and didac- 
tical manners in teaching mathematics. 

As many of the DERIVE Users are also 
using the CAS-77v the DNL tries to com- 
bine the applications of these modem tech- 
nologies. 


Editor: Mag. Josef Bohm 
A-3042 Wiirmla 
D'Lust 1 
Austria 

Phone/FAX: 43-(0)2275/8207 
e-mail: nojo.boehm@pgv.at 


Contributions: 

Please send all contributions to the Editor. 
Non-English speakers are encouraged to 
write their contributions in English to rein- 
force the international touch of the DNL. It 
must be said, though, that non-English 
articles will be warmly welcomed nonethe- 
less. Your contributions will be edited but 
not assessed. By submitting articles the 
author gives his consent for reprinting it in 
the DNL. The more contributions you will 
send, the more lively and richer in contents 
the DERIVE & CAS-77 Newsletter will be. 


Next issue: March 2005 

Deadline 15 February 2005 


Preview: Contributions waiting to be published 

Finite continued fractions St. Welke, GER 

Some simulations of Random Experiments, J. Bohm, AUT & L. Kopp, GER 

Wonderful World of Pedal Curves, J. Bohm 

Another Task for End Examination, J. Lechner, AUT 

Tools for 3D-Problems, P. Luke-Rosendahl, GER 

ANOVA with DERIVE & Tl, M. R. Phillips, USA 

Hill-Encription, J. Bohm 

CAD-Design with DERIVE and the Tl, J. Bohm 

Avoiding Convolution and Transforming Methods, M. Lesmes-Acosta, COL 

Farey Sequences on the Tl, M. Lesmes-Acosta, COL 

Simulating a Graphing Calculator in DERIVE , J. Bohm, AUT 

Henon & Co, J. Bohm 

Pringles, B. Grabinger, GER 

Challenges from Fermat, Bj. Felsager, DEN 

Actuarial Mathematics, M. R. Phillips, USA 

Are all Bodies falling equally fast, J. Lechner, AUT 

Modelling Traffic Density, Th. Himmelbauer, AUT 

Do you know this? Cabri & CAS on PC and Handheld, W. Wegscheider, AUT 
Two Stage Least Squares, M. R. Phillips, USA 

and Setif, FRA; Vermeylen, BEL; Leinbach, USA; Koller, AUT, 

Keunecke, GER and others 

Impressum: 

Medieninhaber: DERIVE User Group, A-3042 Wiirmla, D'Lust 1, AUSTRIA 
Richtung: Fachzeitschrift 
Herausgeber: Mag. Josef Bohm 
Herstellung: Selbstverlag 
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♦♦♦ Describing the many treasures which are more or less hidden among the User 
contributed files. Most of the users don't know about the contents and how to use. 

♦♦♦ It could be useful to deal with one or the other files of the MATH-folder, too. 


At this occasion the editor expressed his thanks to Theresa Shelby, Albert Rich, David Stoutemyer and 
Bernhard Kutzler for their great cooperation. 

The immediate contact between the users and to the responsible people is one of the great advantages 
of Derive and the CAS-TIs and - 1 am sure - a major part of their success. 

Speaking about thanks we must not forget Noor Bohm, who has done all the administration work since 
1991. 

The report closed with an outlook on future contributions. 

Then Bernhard Kutzler reported that he had audited the finances of the DUG and had found that 
everything was ok. 

There was one proposal for the management committee for the next period, which was accepted with 
one voice: 


Josef Bohm 
Barbel Barzel 
Noor Bohm 
Josef Lechner 
Bernhard Kutzler 
Walter Klinger 

Johann Wiesenbauer suggested that due to the fact that the publications are now on the web we could 
pose a "Problem of the Month” to be tackled by our members. So we will put again "Challenge” into 
the next newsletter 64 and wait for the responds. If there is need in more "problems" we can start with 
them at any time, if there is not, it would be lost time and efforts to bring them on the web. 

Lana Moore from TI announced cooperation with TI by promoting the User Group and the Newsletter 
via the TI-Derive-Homepage. Many thanks for this valueable support. 

Thanks also to all friends who came together and honored the meeting with their presence. 

The next meeting was announced for the next DERIVE Conference which will take place in 2006 in 
Dresden, Germany. 


The meeting was closed after 2 hours. 





p 


4 


Derive- and CAS-TI-User Forum 
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Total Differential with DERIVE and the TIs 


Don Pillips 

Thanks for DNL #55! It was great as usual. I really liked how you did my article on TVM comparing 
my routines with yours and the TI-89. 

I've expanded my routines for actuarial math. So, if you're planning on publishing it at some future 
date, please wait for the "new, expanded, greater" version. 

Also, I've noticed, as I'm sure many have, that Derive does not compute total differentials. I've at- 
tached a file which corrects that. Maybe Albert, et.al., will include this functionality in a future version 
of Derive. 

Sincerely, 

Don Phillips 


Example 1: Find the first order total differential for y = x A 2. 

2 

#2: TotalDiffer ent ialfy = a ) 

#3 : dy = 2 -a ■ da 


Example 2: Find the first and second order total differentials of z = 4x A 3y A 2. 

^ 2 

#4: TotalDiffer ent ial(z = 4-a -y ') 

2 2 3 

#5: dz = 12-a -y ■ da + 8-a -y-dy 

5 2 

#6: TotalDiffer ent ialfz = 4-a -y , 2) 

2 2 2 3 2 

#7: d2z = 24-a-y -da + 48 -a -y-da-dy + 8-a ■ dy 


Example 3: Let z = -/(x A 2+y A 2) . Use a total differential to 
approximate the change in z as (x,y) varies from the point (3,4) 
to the point (3.04,3.98). 


2 2 

#8: TotalDifferentialfz = V (a + y )) 

a ■ da y ■ dy 

dz = + 

#9: 22 22 

VO + y ) VO + y ) 


5UE5T 


# 10 : 


a ■ da 


dz = 


y-dy 


2 2 2 2 
VO + y ) VO + y ) 


[a , y , da , dy] , [3, 4, 0.04, -0.02] 


#11: dz = 0.008 


Josef, 

Attached is the total differential program for the TI-89. 

It takes two arguments: the equation and the degree of dif- 
ferentiation. I just wish there was some way to put a default 
value in a TI-89 function or program. 

Regards, Don. 

You can find the TI-89, the TI-92 and the V200 grouped file 
among the files for download. Josef 


Algebra 


FHt Y— FE T Ffir Y '| 

i|Calc|01her|PrgnIQ|ciean Up| | 


' tot_dif(z = 4x 3 y 3 , l) 
dz = 12 ■ dy ■ x 
' tot_dif(z = 4 x 3 y 3 ? 2 ) 
d2z = 24 ■ dy 2 ■ x 3 ■ y + 72 ■ dx ■ dy ■ x 2 ■ y 2 + 24 > 
■ tot_dif(z = , 0 dz = + 


tot-dif <z=JXx A 2+v A 2> ,1>I 
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Note from Josef Lechner and Ove Kroll concerning working with angles in DERIVE 6: 

Setting Mode Degree/Radian has no effect on the input of angles: 

#1 : Angle := Degree 

#2: TAN(45) 

simplified 

#3: TAN(45) 

approximated 

#4: 1.619775190 

but tan(45°) gives simplified and approximated: 

#5: TAN(45°) = 1 

Now the same in Mode Radian: 

#6: Angle Radian 

#7: TAN(45) 

#8: 1.619775190 

#9: TANf45°) = 1 

Setting the Mode has effect for the output of results: 


#10: 

Angle := Degree 




#11: 

SOLVE(TAN(cO = 1, a, 

Real) = i 

= 

225° v a = (-135)° v a : 

#12: 

NSOLVE(TAN(cO = 1, a, 

, Real) = 

c« = 

: -5 ,497787147) 

#13: 

NS0LVE(TAN(O = 1 , ct , Real) : 

= Co 

Lm 

II 

#14: 

Angle := Radian 




#15 : 

Precision := Exact 




#16: 

Notation := Rational 




#17: 

SOLVE(TAN(a) = 1, a, 

Real) = 

r \ 

s 

II 

5 -n 3 -n 

"■ r a — ' f a 

4 4 

#18: 

NSOLVE(TANCa) = 1, a, 

, Real) = 

Co = 

; -5 ,497787147) 

#19: 

NSOLUTIONSfTAN(a) = 

1, a, Real) 


1° 




#20: 

[-315] 




#21: 

ATAN(l) 




#22: 

0.7853981633 




#23: 

ATAN(l) 

= 45 





i° 


But in DERIVE 6 you can use the ARC-functions: 

#24: ARCTANtl') = 45 


See Albert Rich's comment on page 14 



Tania Koller and Students (Illb) Model Reality on the Voyage 200 



See also page 22 
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Quadratic Approximations for Integration 

Giora Mann Nurit Zehavi 

Levinsky College of Education, Israel Weizmann Institute of Science, Israel 

Abstract 

The fundamental theorem of Calculus states that the definite integral over an inter- 
val is the difference of the antiderivative at the two endpoints of the integration in- 
terval. However, we need to know the antiderivative, but in too many interesting 
cases (for example in computing the length of curves) the antiderivative does not 
exist, or we are unable to find it. 

In this paper we present a didactical sequence in which we start by showing that a 
quadratic approximation of the function to be integrated can give a "good” ap- 
proximation of the definite integral. Furthermore, this approximation depends only 
on three values of the integrand (no antidervative is needed). Refining the ap- 
proximation opens the road to Simpson's rule for computing numerically definite 
integrals. 


Introduction 


The motivation for the didactic sequence presented in this paper originated from comments of curious 
high-school students while learning to compute definite integrals: "It is interesting that the definite 
integral for a given function depends only on the values of the primitive function at the endpoints of 
the integration interval; the area under the graph is not affected by the behavior of the primitive func- 
tion and its derivative inside the interval of integration." This intimate connection between the definite 
integral and the given function is expressed in the fundamental theorem of Calculus: 

b 

For F'(x) =fix), | f(x) = F(b) - F(a) . 

a 

The following is a typical problem given to students: Find the area between the graph of the function 
J{x) = 0.5x 3 - 2x 2 + 4 and the x - axis, above it (see Figure 1). We demonstrate the solution and explo- 
ration of the problem using Derive : 



Figure 1 . 
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fGO := 


3LVEC0.5.X - 


x = 1 - J5 v x = J5 + 1 


CO . 5 - x - 2 «x + 4) dx = 


5-45 IB 


l-i/5 


APPROX J (0.5-x - 2-x + 4) dx = 

L 1 - ^5 


5-V5 IB 


= fals* 


Why is Derive' s reaction to the Approx command - "false” ?! 
Let us approximate the left side first and then the right side: 


#6: 8. 060113295 = 


-V5 13 


#7: 8. 060113295 = 8.060113295 


So, why "false”? 

The conflict is resolved when we change the number of digits: 
#8: NotationDigits := 12 


#9: APPROX J (0.5-x - 2-x + 4) dx = 8. 06011329556 

L 1 -75 J 

f 5.^5 18 ^ 

#10: APPROX + = S. O60113295S3 

L 3 3 J 

We see that the software implements different algorithms for approximating integrals and for ap- 
proximating real numbers. It is worthwhile to learn more about the algorithm for approximating inte- 
grals. 

A quadratic approximation for a definite integral 


Let us look at g(x) : = x 3 - 3x + 6 and find the area cir- 
cumscribed by the graph of g(x), y = 0, x = -2, x = 1 


y=x A 3-3-x+6 


Figure 2. 4 -3 j -2 -1 1 


-15 
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We want to fit a quadratic function q(x ) to g(x); 

" T -2 gC-2) T 

#13: qO) := FIT [x, A«x + B «x + c] , -0.3 g(-0.5) 


-0.3 gC-0.5) 

1 gCl) 


3-x 3-x 

#14: qO) := + 7 


Figure 3 indicates that (a) the graphs of the two functions in the integration interval are of different 
types (e.g., only g(x) has an inflection point), and (b) the area under q(x) seems to be the same as the 
area under g(x). 

y / 

y=x A 3-3-x+6 -i r / 


And indeed. 


/ 


/ 


/ 



X 

\ : 


\ 



\ 

7 

\ 


Figure 3. 


i i 

Jg(x) = 18.75 | q(x) = 18.75 


This is true for any definite integral of a polynomial of the third degree. We can show that easily, but 
having CAS at our disposal, we prefer to present a broader framework. It is enough to deal with the 
simplest polynomial of n : 

n 

Pn(x) := A-x 


QnCx) := FIT 


2 

_x , a-x + b-x 


xO - 

k 

A- (xO - k) 



n 

xO 


A-xO 

xO + 

k 

A- (xO + k) 


A(n) := J Qn(x) dx - J Pn(x) dx 

xO - k xO - k 


Computing the difference for n = 1 . . .7 yields: 






plO 


G. Mann & Nurit Zehavi: Quadratic Approximation 


D-N-L#56 


TABLEfAfrO, n, 1, 7) 


4 ■ A ■ k -aO 


4 ■ A ■ k .(2-k + 21-xO 


4 ■ A ■ k -xO-tZ-k + 7-aO 


What conclusions could one make from the table? For n = 1,2 the difference between the integrals of 
the given polynomial and the quadratic approximation is zero, because the approximation is in fact the 
same polynomial. For n = 3, the fact that we get zero difference seems peculiar because we have two 
different integrands. The insight we gain for n = 4 is that the difference depends only on the width of 
the interval, and it is not affected by the location. For n > 4 the difference depends, as expected, on 
both the width and location of the interval. 

Students may wonder why we need to approximate the definite integral of a cubic polynomial by the 
integral of a quadratic function. Well, let us try to find the circumference of the braid created between 
y = sin(x) and y = cos(x) from the first intersection point to the second to the right of the origin (Fig- 
ure 4). 



Figure 4. 


The length of a curve is expressed by the formula \Jl + (y f ) 2 dx . 


So the circumference is twice the integral J ^/l + cos 2 (x)dx , which is not easy to compute without a 


CAS. We can start by approximating such integrals with a quadratic function by computing their coef- 
ficients. But this is one of those moments in mathematics when working generally enables us to see 
surprising results unseen otherwise! 
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Approximating an integral with (only) three values of the integrand 

Assume that q(x) = Ax 2 + Bx + C is a function that meets a given function f[x) at three points: the end- 
points and the midpoint of the interval [a, b\. 

A a ) = q(a), f[b) = q(b) f((a + b)/2) = q((a + b)/ 2) 


f(a) = A-a + B-a + C 
2 

f (b) = A-b + B-b + C 
f 


" a + b ^ 

- 

" a + b ^ 




. 2 j 


. 2 j 



2 

" a + b ^ 

A 

■ a A 




i 


— 

1 

. 2 j 


4 


2 a + b 

+ B + C 


+ C 


b 3 3 2 2 

f 2 A-(b-a) B ■ (b - a ) 

| (A-x + B-x + C) dx = + + C«(b - a) 

a 3 2 

b 2 2 

f 2 (b - a) ■ (2 ■ A- a + 2-A-a-b + 2-A-b + 3-B-a + 3-B-b + 6-C) 

| (A-x + B-x + C) dx = 

a 6 


From the above we conclude that: 


[ q(x) = + 4/(-2— ) + fib)) . 

I 6 2 

The meaning of the result is that the area under the function f(x) over the interval [a, b ] is approxi- 
mately the area of a rectangle whose length is the length of the integration interval and whose height is 
the weighted mean (with weights 1,4, 1) of the three values of the function at a , (a+b)l 2, and b. As 
one can see, the quadratic function “disappeared”. We used only its existence to obtain the approxima- 
tion, which depends only on three values of the given function and the interval of integration. There- 
fore, we reflect again on the fundamental theorem of Calculus that requires the values of the primitive 
function at the endpoints of the integration interval. Clearly, if we are unable to find a primitive func- 
tion, we can easily get a reasonable approximation of the integral, using only the values of the inte- 
grand at the endpoints and the midpoint of interval. 


Toward Simpson’s rule 

We come back to the approximation of the integral needed for finding the circumference of the braid: 

5 n/ 5 Ti / 

/ 4 . i / 4 

J yjl + cos 2 (x)dx . Using the last result for h{x) := -yl + cos 2 (x) , we approximate J h(x)dx by: 


APPROX 


5 ■ 7T 


1 

A 6 

r i r r tt \ 


+ 4 ■ h 


3-7 T 


+ h 


APPROX 


Vo V V J 


3 -77 


A 4 y 
f 3-tt 


5 ■ 7T 


3 . S4 764949 04 3 


+ 4 ■ h + h 

^ 4 J 4 ;;; 


= 1 . 224744S7139 


p!2 
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n \[6 



S-tt/4 


APPROX 


I 


7 t/4 


J(1 + COS 00 ) &a 


Figure 5. 


3 . S201977SS97 


Using Derived numerical integration, we get some control of our result: 

Obviously if we want a better approximation, we should divide the integration interval into sub- 
intervals and approximate the curve by parabolas. This is in fact the idea underlying Simpson's 
method. 

We proceed by dividing the interval into two sub-intervals: 


5 ■ 7t/4 



B -77/4 


5 ■ 7T /4 


J 

h(x) 

dx = 

J 

hOO 

dx + J 

h 00 dx 

7t/4 



7T/4 


B -77/4 


Applying 

: the method, 

we get for each 

sub-interval: 


5 ■ 7t/4 



B ■ 7t/4 


5 ■ 7t/4 


J 

hOO 

dx = 

J 

h(x) 

dx + J 

hOO dx 

7t/4 



7t/4 


B ■ 7t/4 


B ■ 7t/4 



1 r b 

■7T 

7T 'j f f TT 

■\ /■ 


J hO<) dx = ■ 

7t/4 6 


+ 4-h 


J v v h y 


7 r 

v 2 > 


+ h 


3-7 T 






3 -7t/4 

J hO) dx = 1.GSS4724GG2S 

7T/4 


S-tt/4 1 r 5 -it 

J hO<) dx = ■ 

3 ■ 7t/4 6 


3 -7T 


> r r 


J v v 


3 -7T 


+ 4-h(7r) + h 


S ■ 7T 






3 -7t/4 

J hOO dx = 2 . 122235S9443 

3 -7t/4 


Combining the above we get the following approximation: 
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5 - 77/4 

J h&O dx = 1.6SS4724662S + 2 . 12223589443 

77/4 

5 -77/4 

J h(>) dx = 3 . 81070836071 

77/4 

Doubling the number of intervals yields a better approximation: 


5 -77/4 1 f 5 -77 77 

J h(V) dx = ■ - 

77/4 24 l 4 4 J 


4 ■ h 


3-77 


2-h 


4-h 


5 ■ 77 


2-h 


3 ■ 77 


4-h 


7 -77 


s 


+ 2-h(77) + 4-h 


9-77 


5 -77 


5-77/4 

J hOO dx - 3 . 8202S240611 

77/4 

Comparing the three results that we got by applying the method of quadratic approximation with the 
numerical approximation by Derive shows that the greater the number of sub-intervals, the better the 
approximation that we obtain: 

APPR0X( | 3 . 84764949045 - 3 . 82019778897 | ) = 0.027451701471694116 
APPR0X( | 3 . 81070836071 - 3 . 82019778897 | ) = 0 . 0094894282673033272 


APPR0X( | 3 . 82028240611 - 3 . 82019778897 | ) = 8. 46171332SS9279S9 -10 


-5 


Approximating integrals by Derive 

This may be a good time to tell students that Derive actually uses an adaptation of Simpson's rule to 
numerically approximate definite integrals. 


We previously saw that increasing the number of sub-interval increases the accuracy. Simpson's esti- 

mate of j / (x)dx , when the interval is divided into n sub-intervals of equal length h = , is the 

J 2n 


sum: 


h 


[f(x 0 ) + 4f(x 1 ) + 2f(x 2 ) + 4f(x 3 ) + ...2f(x n _ 2 ) + 4f(x n _ l ) + f(x n )], 


where the error is at most — — • m 4 denotes the maximum value of |/ 4 (x)| for x in [a, b]. 

180 

Simpson's estimate is exact for polynomials of the form y = a X* + b x 1 + c x + d, whose fourth deriva- 
tive is zero. The error in using Simpson's methods for other functions involves the fourth derivative. 


We now reflect on the expressions of error fory = Ax , n < 5 that we obtained when we computed the 
error for polynomials. 


For 72 < 4 we got zeros, which agrees with the fact that M 4 is zero. For n = 4 we got 


4 Ak 5 


15 
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1 . . ( b - a)M h 4 

We substitute in : 

180 


M 4 = 4! -A, b - a = 2k, h = k and get: 


2k-4 \-A-k 4 _4 A-k 5 
180 ” 15 


The significance of the last result is that our error for the polynomial y = Ax 4 is exactly the accuracy 
allotted by Simpson's rule. 


Concluding remarks 

The didactical sequence we described above utilizes CAS to make numerical integration by CAS less 
mysterious. We start by a ‘didactic moment’ in which students obtain two different results in using the 
software to approximate definite integrals. Next they realize that for a quadratic approximation, the 
coefficients of the quadratic function are not needed; the approximation is produced by using only 
three values of the integrand. At this stage the road to Simpson’s method is open; moreover when we 
mention the accuracy of the method the CAS can be used to get ‘some feeling’ of the key indicator of 
the error. 


Albert Rich's comment on Angular Mode Settings 

Q: When calling on trig functions, how do I enter angles in degrees? 

A: In Derive 6, the ° operator is used to enter an angle in degrees. The ° operator can be 
entered by clicking on it on the math symbol toolbar, pressing Ctrl+O, or by typing deg on the 
expression entry line. For example, SIN(45°) simplifies to SQRT(2)/2. Unlike earlier ver- 
sions of Derive, selecting Degree in the Angular Unit field of the Simplification tab of the Op- 
tions > Mode Settings command only effects the display of angles, not how angles are en- 
tered. 

Q: In approximate mode, how do I get the inverse trig functions to return angles in degrees 
instead of radians? 

A: In approximate mode, the built-in inverse trig functions (e.g. ASIN, ACOS, ATAN, 
etc.) always return angles in radians, even in Degree mode. For example, in Degree mode 
ATAN(1) simplifies to 45° but approximates to 0.7853981633. To always get angles returned 
in degrees use the inverse trig functions (e.g. ARCSIN, ARCCOS, ARCTAN, etc.) defined in 
MiscellaneousFunctions.mth instead of the built-in functions. For example, ARCTAN(I) sim- 
plifies and approximates to 45. 


Note the new inverse trig functions: ARCSIN, ARCCOS, etc. 
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Applications of the Moore-Penrose Inverse of a Matrix 

Karsten Schmidt, FH Schmalkalden, Germany, kschmidt@fh-sm.de 


Introduction 

After giving an introduction to the Moore-Penrose inverse of a matrix, and its computation in 
DERIVE , in DNL #50 (Schmidt 2003), this paper deals with two important applications of the Moore- 
Penrose inverse. One is a method for solving a system of linear equations, and the other is the compu- 
tation of the Ordinary Least Squares estimator in linear regression. Some familiarity with matrix alge- 
bra as well as basic understanding of the Moore-Penrose inverse of a matrix (as provided in DNL #50) 
is required. 

Computation and Properties of the Moore-Penrose Inverse 

In order to facilitate working with this paper, the definition and DERIVE functions for the computation 
of the Moore-Penrose inverse of a matrix are repeated from DNL #50: 

For any mxn -matrix A there exists a unique matrix with properties related to those of the inverse of a 
nonsingular matrix. This is the Moore-Penrose inverse, denoted by A + , which satisfies the four condi- 
tions (the transpose of A is denoted by A ) 


AA + A = A 

(i) 

A + AA + = A + 

(2) 

[a + A^ =A + A 

(3) 

+ 

II 

+ 

(4) 

Conditions (3) and (4) require both A + A and AA + 

to be symmetric matrices. Note that A + is 


nxm -matrix, i.e. the dimension of A + is equal to the dimension of A . 

The Moore-Penrose inverse of a matrix can be computed in DERIVE with the following two functions: 


MPIV(a) := 

If DIM(a') = 1 

If (a' - a) 4 , 14,1 = 0 
0 • a' 

a’/Ca’ *a)ilil 

"This is not a column vector!" 

MPI(A, APLUS, aj, dt, c, bt, I) := 

Prog 

APLUS := MPIV(A COL [1]) 

3 := 2 
Loop 

If 3 > DIM(A') 

RETURN APLUS 
aj := A COL [J] 
dt := aj 1 - APLUS' -APLUS 

c := (IDENTITY_MATRIX(DIM(A)) - A COL [1, . 3 - 1] -APLUS) -aj 
bt MPIV(c) + (1 - MPIV(c) • c)/(l + dt-aj)-dt 
APLUS := APPEND(APLUS - APLUS -aj-bt, bt) 

3 :+ 1 
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This time A is singular (det(^4) = 0), its inverse A 1 does not exist. Computing the Moore-Penrose 
inverse 


A = 


{ _2_\ 
25 25 

_2_ J_ 

V 25 25 J 


is nevertheless possible and we find that condition (8) is satisfied for system (10): 


AA + b = 


1 2 
2 4 


f i 


25 25 

_ 2 _ A. 

V 25 25 2 




v 10 y 


fl 2^ 
5 5 

2 A 
V 5 52 




v 10 y 




v 10 y 


=6 


As a third example, look at 

^ 1 2^| _ f 5 ^ 

2 4 


A = 

2x2 


; b = 

2x1 


v 15 y 


( 11 ) 


This time we find that condition (8) is not satisfied: 


AA + b = 


( 1 2 \ 
5 5 

f 5 l 


"7" 

2 4 

V 5 5) 

ll5. 


vl4. 


*b 


System (1 1) is therefore inconsistent. 


CHECKSLE(A, b) := 

If A-MPI(A) ■ b = b 
#1: "consistent" 

"NO SOLUTIONS!" 


#2: A := 

#3: b := 

#4: 

#S : A := 

# 6 : 

#7: b := 

# 8 : 

#9: b := 

# 10 : 


1 2 
3 4 J 
S 

10 

1 2 
2 4 


S 

IS 


CHECKSLE(A, b) = consistent 


CHECKSLE(A, b) = consistent 


CHECKSLE(A, b) = NO SOLUTIONS! 


The function CHECKSLE in the above screenshot checks if a system of linear equations is consistent or 
not, and prints the result on the screen. Since there is no unknown clause in the I F-expression, the 
entire (simplified) I F-expression is returned, which can obviously be fairly informative in cases such 
as the last SLE (consisting of matrix A in #5 and vector b in #9). 


If a system of linear equations Ax = b is consistent, its general solution is given by 

f \ 


x-A + b + 


I 




AA 

J 


z 

nx 1 


where z is an arbitrary vector. 


( 12 ) 
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Note that since the vector jeR" in (12) is arbitrary, we might simply choose - () . Consequently, 

one (possibly unique) solution of a consistent system of linear equations Ax = b is always given by 

x = A + b 

The following function SOLVESLE either solves a system of linear equations Ax = b where the matrix 
A and the vector b have been passed as parameters, or displays a message if the system is inconsistent. 


z := VECTOR(VECTOR(APPEND(z, J), i, 1), 3, 1, DIM(A’)) 

SOLVESLECA, b) 

If A • MPI(A) . b = b 

MPI(A) . b + (IDENTITY JV|ATRIX(DIMENSION(A')) - MPI(A) . A) . z 
"No solution(s) ! " 


We now want to compute the solution(s) of the above three systems. We start with system (9). Since A 
is a nonsingular matrix, we have A + = A 1 , and (12) simplifies to 

x - A + b + (j - A + Al^z 


A x b + 


= A~ l b 


I-A 'A \z 


for any choice of z- Obviously, the general solution (12) simplifies to a unique solution if A is 
nonsingular. Hence 


x = A l b 


r -2 

l ^ 

(5) 


f°l 

3 

1 

yov 

— 

5 

V 2 

2y 


\2j 


For system (10) we get 


x = A 


r b + (l-A + A)z ■■ 




V 5 


jz 1 +jz 2 + 2 


In this case there is an infinite number of solutions. For example, by choosing z = 0 we get 


x = Ab = 


v^y 


#1: A :z 


1 2 
3 4 
5 

10 


SOLVESLECA, b) = 


1 2 
2 4 


SOLVESLECA, b) = 


4-zl 2- z2 


S S 

2-zl z2 


+ 2 


# 8 : 


15 
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The above screenshot shows the capability of the function SOLVESLE to handle all three possible sce- 
narios in considering a system of linear equations: a unique solution, an infinite number of solutions, 
and the case that no solution exists. 


Linear Regression and the Moore-Penrose Inverse 

We consider the (multiple) linear regression model 
y = X p + u (13) 

Nx 1 NxKKx 1 Nx 1 

where y is the vector of observations on the dependent variable, X the regressor matrix, /? a vector of 
parameters, and u a vector of disturbances. 

Denoting an estimator of the unknown parameter vector /? by (5 , we have 

y = Xp 

u = y-y 

where y is the estimate of y using ft , and u is the vector of residuals. 

The most popular estimator for / 3 is the (Ordinary) Least Squares estimator which minimizes the sum of 
squared residuals 


?>(/?) 




Note that 


<p(p) = (y-Xfi)(y-Xfi) 

= y'y - y'Xp - p'X'y + fiXXfi 
= P'XXp-lyXfi + y'y 

is a convex function since XX is a nonnegative definite matrix. Therefore, finding its first derivative 

= p'[xx + (xx)\- 2 yX 


= 2pXX-2yX 

and setting it equal to 0 is necessary and sufficient to determine the minimum of : 
2B , XX-2y , X= 0 <=> XXB-X*y= 0 XXB = X'y 

1 xK Kx 1 

The last equation constitutes the so-called system of normal equations. 

Under the (usual) assumption that rank(X) = K , which assures that XX is nonsingular, we can eas- 
ily derive the Least Squares estimator from the normal equations 

XXj3 = X'y <=> (XXy 1 XX fi = (XXy l X'y <=> j3 = (XXy l X'y 

i 

One might think that the system of normal equations is inconsistent if rank(X) < K . However, this is 


not true. 
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failures 

t [temp, in °C] 

temp, in °F 

failures 

t [temp, in °C] 

temp, in °F 

3 

12 

53 

0 

21 

70 

1 

14 

57 

0 

21 

70 

1 

14 

58 

0 

22 

72 

1 

17 

63 

0 

23 

73 

0 

19 

66 

2 

24 

75 

0 

19 

67 

0 

24 

75 

0 

19 

67 

0 

24 

76 

0 

19 

67 

0 

24 

76 

0 

20 

68 

0 

26 

78 

0 

21 

69 

0 

26 

79 

1 

21 

70 

0 

27 

80 

1 

21 

70 

0 

27 

81 


The following screenshot shows an algebra window and a 2D-plot window. The first three expressions 
in the algebra window are the contents of the file T4Data.mth. Expression #1 defines a 24x2 -matrix 
Data, which was entered in transposed form to save space (unfortunately, this requires simplification 
of Data prior to plotting the points in the 2D-plot window). In expressions #2 and #3 the data are rear- 
ranged according to the definition of the linear regression model (13). X denotes the regressor matrix, 
containing a column of ones (for the y-intercept), and a column with the observations on the independ- 
ent variable (temperature in Celsius), y denotes the vector of observations on the dependent variable 
(number of failures). 

Expression #4 is the formula for the computation of the Least Squares estimator using the Moore- 
Penrose inverse. Approximating #4 yields #5, which is in turn used to define and finally plot the 
straight line which is the result of the Least Squares estimation. 


xj 


File Edit Insert Set Options Window Help 


^injxj 



| Algebra 2 T4Data.mth 




# 1 : 


Data := 


12 14 14 17 15 IS 15 15 20 21 21 21 21 


3111000000 
21 22 23 24 24 24 24 26 26 27 27 

00020000000 
#2: X := APPENDJZ0LUMNS(VECT0R (VECTOR (1, i, 1), i, DIM (Data)) , Data 

COL [1]) 

#3: y := [Data COL 2]’ 

#4: MPI(X)-y 

2.628625402 
0.1051227914 


EB Cross: -0.9677419 , 2.733871 |Center: 10 , 1 .25 

[Scale: 5 : 0.5 

jj „ = m * * x r 


II ildll il ddddddddddd d UU4 

Ililxj 
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Considering the relatively high coefficient of determination ( R 2 =0.3), and the fact that the slope 
parameter is statistically significant (at the 1%-level; both values not shown in the screenshot), the 
above result is fairly reliable. 

Since the pre-launch ambient temperature on January 28, 1986, was -1 °C (31°F), the prediction from 
the above regression would have been 

failures = 2.629 - 0. 105 (-1) = 2.734 

i.e. 2 or 3 O-ring failures were to be expected according to our regression result. Nevertheless, the 
space shuttle was launched. Less than two minutes into the flight, due to O-ring failure, leaking fuel 
was ignited by a rocket engine, and Challenger exploded. 


Reference 

Schmidt, K. (2003), An Introduction to the Moore-Penrose Inverse of a Matrix, The DERIVE- 
Newsletter #50 (June 2003), 12- 18. 


Artif icial and Real Fountains 




Picture: Fountain in St. Petersburg, Russia (Tania Koller) 
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Kvadratisk programmering 

Et eksempel pa et valgfrit emne i symbolsk matematik 

Quadratic Programming 

An Example for an Optional Topic in Symbolic Mathematics 

Bjoern Felsager, Denmark 


Netop fordi symbolske programmer arbejder med vilkarlige udtryk kan man lige sa nemt lave kvadra- 
tisk programmering, som lineaer programmering (eller en hvilken som heist anden form for program- 
mering!). Det er jo alligevel de samme faciliteter, man traekker pa! Her vil vi se pa nogle eksempler pa 
kvadratisk programmering hentet fra en standardlaerebog for handelsgymnasiet (Soren Antonius et 
al.): 

En produktion af to varer A ogB er underlagt betingelser, som kan udtrykkes ved folgende 
uligheder, hvor x er antal enheder for A ogy er antal enheder for B: 

Production of goods A and B underlies some restrictions which can be expressed by 
the following inequalities for number of units x of A and number of units y for B. 

2x + 3 y < 240 
2x + 2 y < 180 
4x +y < 240 
x > 0 

y>0 

Af de tre Ibrste betingelser folger, at de kritiske x-vaerdier (dvs. skaeringen med x-aksen) er 120, 90 og 
60, ligesom de kritiske y-vaerdier er 80, 90 og 240. Altsa skaerer polygonomradet x-aksen i 60 og y- 
aksen i 80. Vi starter derfor med at vaelge graffummet 


-10<x<70 og -10<y<90 . 

Vi indskriver derefter kriterieme 


2-x + 3-y < 240 a 2-x + 2-y < ISO a 4 ■ x + y < 240 AX>0Ay>0 


og tegner kriterieomradet: 


The restrictions describe a region which has 
its boundaries in critical x- and y-values. 

The intersection points of the boundary-lines 
result in a polygon. 

The points on the axes are easy to find. But 
we need also the intersection points of the 
restriction lines. 
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[rl := 2-x + 3-y = 240, r2 := 2-x + 2-y = ISO, r3 := 4-x + y = 240] 

SOLUTIONS (rl A r2, [x, y]) = [[30, 50]] 

SOLUTIONSCrl a r3, [x, y]) = [[4S, 48]] 

SOLUTIONS (r2 a r3 , [x, y]) = [[50, 40]] 

Det kan ogsa betale sig at indskrive ligningeme 
for radomradet, sa vi fx bestemme skaerings- 
punkteme: 

Ved at benytte Solutions i stedet for Solve far 
vi skaeringspunkteme ud pa koordinatform, sa de 
kan tegnes direkte: 

By using Solutions instead of Solve we ob- 
tain the intersection points in coordinate form 
and the points are immediately ready for plot- 
ting. 

Vi har styr pa kriterieomradet! 


Under passende antagelser bliver den samlede omscetning nu givet ved det folgende kvadratiske udtryk 
ixogy: 

/(x,y) = 20x-y + 25y-^p 

Det er altsa denne funktion vi skal maksimere i polygonomradet, sa vi ser pa nogle niveaukurver , fx 
niveaukurven gennem (25,25), dvs. kurven med ligningen / (x, y) =f (25, 25) . Faktisk kan vi lige sa 
godt tegne en familie af niveaukurver / (x, y) = k , sa pa basis af vaerdien af / (25,25) gaetter vi pa et 
passende interval af familieparameteren k. 



Under certain circumstances revenue is 
given by the following quadratic expression: 

fix, y) = 20x - y + 25y - 

We try to find the maximum value for this 
function without violating the restrictions. We 
inspect so called level curves (contour curves 
of the surface f), eg the level curve contain- 
ing point (25,25). This is a curve with equa- 
tion f(x,y ) = f (25,25). But we can also create 
a whole family of level curves f ( x,y ) = k with 
appropriate interval for parameter k. 



Now in times of slider bars we can present the family of level curves 
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12 12 

#6: f(x, y) := 20 - x x + 25 -y y 

3 4 


#7: f (25 , 25) = 760.4166666 

Plot the next curve in blue 


12 12 

its: 20 -x x + 25 -y y = f (25 , 25) 

=: 4 


Plot the family of curves in grey 

#9: VECTOR (f(x, y) = k , k, 500, 1000, 25) 

#10: f(x, y) = k 



Det kunne godt se ud som om niveaukurverne er ellipser, samt at centrum for disse ellipser ligger inde 
i kriterieomradet. Vi kan bakke analysen op med et tredimensionalt billede af grafen for omsaetnings- 
funktionen sammen med niveaufladen for (x, y) = (25,25), dvs. z = 760,41666. . . 

The level curves are ellipses with their center within the feasible region. We can perform a 
3D-analysis plotting the revenue surface together with the level plane for (x,y) = (25,25), 
i.e. z = 760.416... (and we can also introduce a slider bar for level planes z = k). 



Det er svaerere at fa kriterieomradet med, fordi kriterieligningerne giver anledning til lodrette planer. 
Vi ma derfor indskrive dem pa parameterform (og saette parameterintervalleme fornuftigt!): 

We want to add the boundaries which appear as vertical planes. For plotting them we have 
to rewrite them in parameter form and set meaningful intervals for the parameters. 


Planen med ligningen 2x + 3 y = 240 tegnes altsa med parameterfremstillingen 

2(120-x) 


x. 


3 


z 


0 <x<70, 0 <z < 1000 
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SOLVE (rl, y) = 



2-C120 


3 


x) 






SOLVE (r2, y) = (y = 90 - x) 


SOLVE- r3, y) = Cy = 240 - 4-x) 


2 -(120 - x) 

x, 

3 


[x, 90 - x, z] 


[x s 240 - 4-x, z] 



Vi kan beregne toppunktet for det kvadratiske polynomium / (x, y) ved at lose toppunktsligningerne i x 
henholdsvis y. Det gores nok nemmest ved at differentiere: 

DIF ( f (x , y ) , x) =0 og DIF ( f (x , y ) , y ) = 0 

Ved at tese disse to ligninger med hensyn til x og y finder vi altsa den optimale produktion: 


We can find the top of the surface f ( x,y ) by means of calculus and solving the respective 
equations we obtain the optimal production plan. 


■VJ 


d 


SOLUTIONS 


— f(x, y) = 0 a — tVx, y) = 0, 
vdx dy 


[x, y] 


[[SO, 50]] 


f (30 , 50) = 925 
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Som forventet ligger toppunktet inde i kriterieom- 
radet og den maksimale omssetning er altsa givet 
ved /( 30,50) = 925 . 

The top of the surface lies within the feasible 
region and the maximum revenue of 925 will 
be reached by producing x = 30 units of A 
and y = 50 units of B. 


Bemcerkning: Vi kan ogsa finde toppunktet mere geometrisk ved at se pa skaeringen mellem niveau- 
kurverne og en familie af vandrette henholdsvis lodrette linjer. Vi loser altsa ligningen for niveau- 
kurveme / (x, y) = k med hensyn til y (henholdsvis x): 

Comment: One can find the top point by geometric means only. We solve the general level 
curve f(x,y) = k for y (and for x): 


S0LUTI0NS(fO, y) = k, y) 


2 2 
2-/3-C/C- x + 60 -x - 3 - (k - 625)) + 25/3) 2-/3 - C25 -^3 - ,/(- x + 60-x - 3-(k - 625))) 


2 2 
2--/3-./C- x + 60-x - 3-(k - 625)) 2 - /3 - /C- x + 60-x - 3 ■ (k - 625)) 

+ 50, 50 


SOLUTIONSfffx, y) = k, x) 


2 2 
/3 - C/C- y + 100-y - 4-Ck - 300)) + 20/3) ./3-C20/3 - ./(- y + 100 -y - 4-Ck - 300))) 


2 2 
^3 /(- y + 100-y - 4-Ck - 300)) V^/C- y + 100-y - 4-Ck - 300)) 

+ 30, 30 


Efter en passende expand af lnsningsudtrykket ses da netop at y-vaerdieme ligger symmetrisk omkring 
50 og tilsvarende fas at x-vserdierne ligger symmetrisk omkring 30. Tegner vi linjerne x = 30 ogy = 50 
ser vi da ogsa netop, at der er tale om symmetriakseme for ellipseme: 

Applying the appropriate EXPAND-command on the solutions we recognize symmetry of the 
curves with respect to x = 30 and y = 50. 
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There is an additional way to visualize the region of possible solutions: 

IF(2-x + B-y < 240 a 2-x + 2-y < ISO a 4-x + y < 240 a x > 0 a y > 0, f (x , y) , 0) 



Dette afslutter det Ibrste eksempel, som var simpelt, netop fordi toppunktet la inde i kriterieomradet. 
Vi folger op med et nyt eksempel, hvor toppunktet ligger udenfor kriterieomradet! 


The result of our first example was not so difficult to find because the vertex of the surface is 
within the boundaries. We proceed with an example with the vertex lying outside of the criti- 
cal region. We just change the revenue function on one place. 


Eksempel 2: Vi prover sa at se pa, hvad der sker, hvis vi sendrer omssetningsfunktionen en anelse til 

f(x,y) = 20x-y + 25y-^j- 

Det aendrer niveaukurveme, sa de nu har centrum udenfor polygonomradet: 
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r d 


d 


SOLUTIONS 


— fl(x, y) = 0 a — tl(x, y) 
^dx dy 


0, [x, y] 


[[ 50 , 50 ]] 

The vertex is outside of the boundaries. 


I stedet skal vi derfor bestemme maksimumspunktet, som det punkt pa randen af polygonomradet, 
hvor niveaukurven netop tangerer polygonomradet, dvs. hvor niveaukurven slipper kriterieomradet. 
Forst ma vi da finde ud af hvilke af begrsensningsfiinktioneme der er tale om. En trace viser, at det er 
den anden af kriteriefunktioneme (2x + 2y = 180): 


For finding the maximum point we have to 
find a point on the borderlines of the feasible 
region (- the circumference of the polygon) 
which is also a point of an osculating level 
curve. At first we can try to find an approxi- 
mative solution by using Trace mode. We 
trace along the second constraint line and 
estimate the osculating point of a possible 
level curve. 

One estimate could be (44.3, 45.7). 



We proceed by applying implicit differentiating function fl (x,y) which gives the slope in any 
point. We try to find a level curve which has not only a point with 2x + 2y = 180 in common, 
but also the slope in this point. 



Randkurven har altsa ligningen y = 90 - x, og dermed haeldningen -1. For at finde hseldningen af ni- 
veaukurven differentierer vi ligningen for niveaukurven, idet vi opfatter y som en funktion af x. 
Hseldningen for kurven er altsa givet ved kommandoen: 

!MP_DlF(fl(x , y) , x , y) (with x,y as default settings) 


Vi finder da: 
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4-Cx - 50) 

IMP_DIF(fl(x, y)) = 

5 ■ (50 - y) 


SOLUTIONS Cr2 a IMP_DIF(flCx, y)) = -1, [x, y]) 


T 400 

410 1" 

j 

_L 9 

9 J_ 


fl 


400 


410 ) 


9 J 


10025 


9 


fl 


400 


410 


9 J 


1113 . SSS38S 


Dermed har vi fundet en oplagt kandidat til den produktion, der giver den maksimale omsaetning, som 
altsa viser sig at vaere: 

x = 400/9 og y = 410/9, 
med den maksimale omsaetning givet ved 

/(400/9, 410/9) = 1113.8888.... 

Vi kan checke losnmgen grafisk ved dels at tegne det optimale punkt, dels tegne den tilhorende ni- 
veaukurve: 


So we found a candidate for the maximum revenue in x = 400/9 and y = 410/9 with gaining a 
revenue of 1 1 13.9. We can check this easily by plotting the respective level curve. 


10025 

flCx, y) = 

9 



70 
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0velser til kvadratisk programmering 

Bemcerkning : En klassisk introduktion til kvadratisk programmering kan man fx finde i 

Mogens Ditlev Hansen : Matematik - Okonomi - Optimering (Abacus 1987) 

Der er gode diskussioner/eksempler i kapitel 4 samt en del supplerende opgaver, fx 

Ovelse 1: Omsaetningsfunktionen (kriteriefunktionen) for en virksomhed er givet ved 

/ (x, y) = -3x 2 + 1 8x - 2 y 2 + 28 y 

hvor produktionen er underlagt betingelserne: 

2x + 5 y < 45 
5x + 2y< 60 
x >0 ,y >0 

Tegn produktionsomradet samt en familie af niveaukurver for kriteriefunktionen. 

Bestem maksimum for kriteriefunktionen. 

Ovelse 2: Omsaetningsfunktionen (kriteriefunktionen) for en virksomhed er givet ved 

/ (x, y) = -x 2 + 1 8x — 3 y 2 +36 y 

hvor produktionen er underlagt betingelserne: 

x + 3y < 21 
5x + 3 ; < 35 
x >0 ,y >0 

Tegn produktionsomradet samt en familie af niveaukurver for kriteriefunktionen. 

Bestem maksimum for kriteriefunktionen. 

Ovelse 3: Omsaetningsfunktionen (kriteriefunktionen) for en virksomhed er givet ved 

f(x,y) = -\ Ox 2 + 1 20x -\0y 2 + 220 y 

hvor produktionen er underlagt betingelserne: 

x + 2 y < 18 
2 x + y < 18 
x >0 ,y>0 

Tegn produktionsomradet samt en familie af niveaukurver for kriteriefunktionen. 

Bestem maksimum for kriteriefunktionen. 

I like this contribution very much and would like to add two more worked examples (Josef). 
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Example 1: 


f{x, y) = 3x + 4y - x 2 + 2 xy - 2 y 2 = Maximum 
x + 2y<7 

2y-x<4 

x,y> 0 


x + 2-y < 7 a -x + 2-y < 4 a x > 0 a y > 0 

2 2 

f3(x, y) := 3-x + 4-y - x + 2-x-y - 2-y 
VECTOR (f 3 (x , y) = k, k, 0, 20, 1) 



S0LUTI0NS(x + 2-y - 7 a IMP_DIF(f3(x , y) , x, y, 1) = IMP_DIF(x + 2-y - 7 , x , y , 1) , [x, y]) 

[[3, 2]] 

2 2 

3-x + 4-y - x + 2-x-y - 2-y = f3(3, 2) 


f3(3 , 2) = 12 



Solution curve in blue 


D 
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7 - x 

X, , z 

2 

4 + x 

x, , z 

2 

z = f3(3. 2) 


IF(x + 2*y < 7 a -x + 2* y < 4 a 
x > 0 a y > 0, f3(x, y), 0> 



Maximum is 1 2 when for x = 3 and y = 2. 

Do you miss a Minimumproblem? Here it is: 


Example 2: 


2 2 

X V 

f(x,y) = - + ^--x-2y = Minimum 

2x + 3y < 6 
x + 4y<5 
x,y> 0 


2 2 

x y 

f 4 (x , y) := + - x - 2 ■ y 

2 2 

2«x + B«y < 6 a x + 4«y <5AX>0Ay>0 
For first information plot two level curves! 
f4(x, y) = o 
f4(x, y) = -1 



2.5 3 3.5 
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We make the 2 nd boundary line to a tangent of one circle of the family of circles with centers 
in (1,2) which form the level curves. We could do this without calculus, too. 


S0LUTI0NS(x + 4-y - S a IMP_DIF(f4(x , y), x, y, 1) = IMP_DIF(x + 4-y - B , x , y , 1) , [*, y]) 


T 13 

IS 1" 

i — 
i — 

17 J. 


f4(x, y) = f4 


r 13 


L 17 


IS '' 
17 , 




Minimum = -2.03 for x = 0.765 and y = 1 .059. 

Reference: The Operational Research Problem Solver, REA, New York 1985 
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A Tool for Generating Tree Diagrams 

Lorenz Kopp, Neumarkt, Germany 


Base figure is a circle /o' with radius rand center M_. 


r := 0.2 


2 2 2 
k r (m_) := (x - m_ ) + (y - m_ ) = r 

1 2 


Meaning of the variables: [ax, ay] starting point, dx, dy increase of x and y to the next circle 
(1 st point right on the top = slope triangle). Radius r is considered. 


Two Branches: one, two, three or four experiments 


tree21(ax, ay, dx , dy) 


ax + r ay 

ax + dx - r ay + dy 
ax + r ay 

_ L ax + dx - r ay - dy 


kr([ax+dx, ay+dy]) 


kr([ax+dx, ay-dy]) 


tree22(ax, ay, dx , dy) 


tree23(ax, ay, dx , dy) 


tree24(ax, ay, dx , dy) 


tree21(ax, ay, dx , dy) 


tree21 


ax + dx , ay + dy , dx , 


tree21 ax + dx , ay - dy, dx , 
tree21(ax, ay, dx , dy) 


tree22 


ax + dx , ay + dy , dx , 


tree22 ax + dx , ay - dy, dx , 
tree21(ax, ay, dx , dy) 


tree23 


tree23 


ax + dx , ay + dy , dx , 
ax + dx , ay - dy , dx , 


dy ^ 

2 j 

dy ^ 

2 j 

dy ^ 

2 j 

dy ^ 

2 j 

dy ^ 

2 j 

dy ^ 

2 j 


The first two generations: 

tree21(0, 0, 2, 2) 



tree22(0, 0, 2, 2) 
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Some hints for plotting the trees: 

Window Tile Vertically, Plot region 0 < x < 8, -4 < y < 4, Set Cross 4,0 and Cross on Center. 
Switch off Axes, Labels, Grids. 

My tip: Set Option Display Grids Intervals 14,8 (Josef) 

See tree23(0,0,2,2) with comments (annotations): 
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We provide some other fundamental types of trees: 


Two generations with three branches. 


tree32(0 , 0 , 2 . 5 , 2 . 5) 



1 st generation 3 branches, 2 nd generation 
2 branches: tree3121(ax,ay,dx,dy) 

tree3121(0 ,0,3,3) 



Finally you can use this tool to create your 
special tree. 


Assume you would like to have a decision 
tree according to the sketch! 


1 st generation 2 branches, 2 nd generation 
3 branches: tree2131(ax,ay,dx,dy) 

tree2131(0 , 0 ,1. 5 , 2 . 5) 
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Warning: More solutions may exist 

I received an email concerning the various different answers when solving equations using different 

CAS-calculators, different OS, different modes. In many cases you find a warning: More So you 

never can be sure, if there are other solutions in addition the shown ones. 

Some time ago I published a DERIVE-tool to overcome this problem for many cases. So I adapted this 
program for the handheld devices and made of MYSOLUTIONS a function mysol(equation, variable, 
lower_bound, upper_bound, step). The equation which caused the discussion was e x = x 4 . 

For my experiments I took another equation which has much more solutions. 

Compare the sol ve-output with the mysol -output: 


23 


Rlgebra|Calc|Qther|FrgmIQ|Clean 


i 


■ nysolU x = x 4 , x, "10, 10, l) 

£ -.315553 1.429612 3.613169} 

■ solueU. 5 ■ sin(3 ■ x) = 2 ■ " ' 3 ' x + 1 , x) 

x = 34.743304 or x = 31.490636 or x = 3 . (► 

Horning Mort solutions may txist 





“ r : ' " t ft y rrr 
; c ■; - |Pr gm 1 0 |C L =;■■= =; •= 



■ nysol(f x = x 4 , x , "10, 10, l) 

£-.315553 1.429612 3.613169} 

■ solue(4. 5 ■ sin(3 ■ x) = 2 ■ £■ " ' 3 ' x + 1 , x) 
^30204 or x= "1.423023 or x = "1.670290 


■ nysol(4. 5 ■ sin(3 ■ x) = 2-f -,3 ' x + l,x, "2> 


mvsol£4-5ttsin<3ttx>=2tte A <~ -3ttx... 

MhIH F;hD HUTD FUNC 1/ 5 


Josef 
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Titbits from Algebra and Number Theory (29) 

by Johann Wiesenbauer, Vienna 


In the first place, as for the last column of this series, which dealt with Josephus per- 
mutations from a programming point of view, I want to apologize to Stefan Welke again 
for having overlooked that he had already written a treatise on exactly the same topic. 
In fact, using a new feature of delete( ) in Derive 6, he could achieve an impressive per- 
formance, and there was little to add as you can read yourself in the accompanying De- 
rive file to the Titbits(28). 

Nevertheless, I hope that the derivation of the nice formula for the "survivor" in the 
special case s=2 depending on n, was interesting enough. By the way, I doubt, if there is 
a similarly simple formula for the general case, but I'm convinced that there should be a 
way to write a very fast program to compute the number of the survivor for general 
numbers n and s. In view of the fact that I have not yet programmed it myself, it may 
be a little audacious, but nevertheless I would like to pose this as a challenge for those 
intrepid readers, who "fear neither death nor devil", when it comes to programming in 
Derive. 

Speaking of a challenge, this reminds me of the nice problem posed by Steven Schone- 
feld in the ACDC column in the DNL #55, p40. As you can read there, this problem 
served as an introductory challenge for people who wanted to get hired by Google as 
engineers. In fact, as you can easily check yourself this problem has been discussed on 
many newsgroups and many solutions in different languages were presented on this oc- 
casion. What follows is my solution of this nice problem in Derive 6. 

PHmeSearch(c , d := 10, n := 1, s := 1000, k_ := 0, 1_, p_, u_ := []) := 

Prog 

s := APPROX(STRING(c - FLOOR (c)) , s + d - 1) 
s := REST (REST (s)) 

1_ := DIM(s) - d 
Loop 

If k_ > 1_ 

RETURN REVERSE(u_) 

p_ := C0DES_T0_NAME(NAME_T0_C0DES(s 4. [1 d])) 

If PRIME?(p_) a DIM(p_) = d 
Prog 

u := ADJ0IN([k_ + 1, p_] , u_) 

n 1 
If n = 0 

RETURN REVERSE (u_) 
s := REST (s) 
k_ :+ 1 
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And here a few remarks as to the programming: 

Basically, the program searches for the first n occurrences of d-digit primes in the 
decimal representation of the positive real number c, which is assumed to be irrational, 
where only the first s (by default 1000) decimal places after the decimal point are used. 
Here the primes should have really d digits, i.e. leading zeros are not allowed. Further- 
more, blocks may exceed that bound of s digits as long as their start is within this 
bound. As a nice consequence, always s blocks are checked for primality not depending 
on the value of d and leading zeros are not allowed, hopefully all in perfect agreement to 
what Steven demanded in his challenge. Ok, maybe not quite, because I set n:=l by de- 
fault, which yields only the first occurrence of a d-digit prime, while Steven seemed to 
be more interested in all occurrences within the given bound, which you get by setting 
n:= inf. Furthermore, the routine returns also the locations of the starting digit of the 
primes, but it's easy to suppress this information, if you are not interested in it, as 
shown in the examples below. 

VECTOR(TABLE(PHmeSearch(c, d) , d, 1, 10), c [e, rr, ^2]) 

T 1 [[1, 7]] II" 1 [[4 ' 5]] II" 1 CC4, 2]] T 

2 [[1, 71]] 2 [[2, 41]] 2 [[1, 41] ] 

3 [[4, 2£L]] 3 [[7, 653]] 3 [[3, 421]] 

4 [[14, 4523]] 4 [[2, 4159]] 4 [[7, 5623]] 

5 [[24, 74713]] 5 [[1, 14159]] 5 [[7, 56237]] 

6 [[12, 904523]] ' 6 [[9, 358979]] ' 6 [[5, 135623]] 

7 [[20, 6028747]] 7 [[3, 1592653]] 7 [[6, 3562373]] 

8 [[64, 72407663]] 8 [[33, 28841971]] 8 [[3, 42135623]] 

9 [[19, 360287471]] 9 [[29, 795028841]] 9 [[4, 213562373]] 

_L 10 [[99, 7427466391]] J L 10 [[4, 5926535897]] J L 10 [[1, 4142135623]] J. 


PrimeSearch(e 

, 10, oo) 1 








99 

123 

149 

171 

182 

201 

214 

218 

254 

. 7427466391 

7413596629 

6059563073 

3490763233 

2988075319 

1573834187 

7021540891 

5408914993 

64S0016S47 

295 

309 

313 

322 

399 

438 

440 

473 

505 

9920695517 

1838606261 

6062613313 

3845830007 

1692836819 

4425056953 

2505695369 

5490598793 

1782154249 

507 

516 

536 

537 

547 

577 

599 

617 

645 

8215424999 

9229576351 

9519366803 

5193668033 

1S252SS693 

8294 8879 3 3 

1730123819 

4039701983 

4804295311 

667 

669 

687 

702 

705 

709 

728 

743 

781 

8194558153 

9455815301 

1332069811 

6181881593 

1881593041 

5930416903 

1934580727 

385 8942287 

4841984443 

819 

856 

870 

892 

901 

921 

- 



1978623209 

3140934317 

3640546253 

8887 07 016 7 

7683964243 

4563549061 
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(PnmeSearch(e , 10, «>)) COL 1 

[99, 123, 149, 171, 182, 201, 214, 218, 254, 295, 309, 313, 322, 399, 438, 440, 473, 505, 507, 516, 536, 537, 547, 
577, 599, 617, 645, 667, 669, 687, 702, 705, 709, 728, 743, 781, 819, 856, 870, 892, 901, 921] 


(Pnme5earch(e , 10, ®)) COL 2 

[7427466391, 7413596629, 6059563073, 3490763233, 2988075319, 1573834187, 7021540891, 5408914993, 6480016847, 
9920695517, 1838606261, 6062613313, 3845830007, 1692836819, 4425056953, 2505695369, 5490598793, 1782154249, 

8215424999, 9229576351, 9519366803, 5193668033, 1S252SS693, S294SS7933, 1730123819, 4039701983, 4804295311, 

8194558153, 9455815301, 1332069811, 61S1SS1593, 1881593041, 5930416903, 1934580727, 3S5S9422S7, 4841984443, 

1978623209, 3140934317, 3640546253, 8887070167, 7683964243, 4563549061] 


DIM(PrimeSearch(e , 10, ™)) = 42 


In particular, these computations show that there are 42 10-digit primes which start 
not later than the lOOO-th digit after the decimal point. Is this in agreement with what 
we had expected? Sure! After all, due to the prime number theorem the density of 
primes near a positive real number x should be about l/ln(x). Hence, we would expect 
about 1000/(10 In 10) « 43.43 10-digit primes for s=1000, which is not too far away 
from the actual number above, indeed! More generally, the dependency of the number 
of d-digit primes on s and d can be roughly described by the function 


1 s 
LN(10) d 


Closely related to the problem above is the problem of finding the starting location(s) 
of any string d of digits within a given positive irrational number c, e.g. e, n , V2 , etc. In 
fact, it is not very difficult to modify the program above accordingly. 

StHngSearch(c , d, n := 1, s := 1000, d_, k_ := 0, 1_, p_, u_ := []) := 

Prog 

d_ := DIM(d) 

p_ := APPROX(STRING(c - FLOOR(O), s + d_ - 1) 

s := REST (REST (ITERATE (APPEND(x_, n 0 n ), x_, p_, s + d_ + 1 - DIM(p_)))) 
1_ := DIM(s) - d_ 

Loop 

If k_ > 1_ 

RETURN REVERSE (u_) 

p_ := sj,[l d_] 

If p_ = d 
Prog 

u_ := ADJOIN (k_ + 1, u_) 
n 1 
If n = 0 

RETURN REVERSE (u_) 
s := REST (s) 
k_ :+ 1 
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There is a small subtlety, you should be aware of, though: Trailing zeros, which had been 
cut off by Derive when approximating c-floor(c), must be appended again. This is done 
by the line 

s := REST (REST (ITERATE(APPEND(x_, "0"), x_, p_, s + d_ + 1 - DIM(p_)))) 

in the program above. (Note that this cutting off of trailing zeros didn't bother us in 
the first program, since we were looking for strings representing primes only!) 

Well, what about the distribution of 0-9 within say the first 10000 digits of n after 
the decimal point? Here you are! 

TABLE(DIM(Str-ingSearchO, STRING(d) , 10000)), d, 0, 9)’ 

'0 1 23 4 5 6 7 S 9" 

. 968 1026 1021 974 1012 1046 1021 970 94S 1014 . 

Now that the year 2004 is coming to an end (hopefully it was a good year for you and 
your family!), we could ask whether it is contained in the decimal representation of n 
starting not later than 10000 digits after the decimal point? Wanna bet? 

Stri ngSearch(pi , "2004", 1, 10000) = [7235] 

Yes! By the way, this bet hasn't been extremely risky, as the chances of a failure are 
only about 1/e « 36.79%, assuming that n is normal, i.e. its decimal places form a 
"good" pseudorandom sequence. (Most mathematicians believe this to be true, although 
no rigorous proof is known!) 

Let's turn to a completely different topic now, which should have already been part of 
my talk in Montreal, but it turned out that there was not enough time for it. It deals 
with the so-called discrete logarithm problem (DLP for short) with many important ap- 
plications in modern cryptography (Diffie-Hellman key exchange, ELGamal cryptosys- 
tem, DSA etc). 

In its most general form it can be stated as follows. Given a finite cyclic group & of or- 
der n, a generator a of G and an element p eG, find the unique integer x, 0 < x < n, such 
that a x =p. This integer x is called the discrete logarithm of p to the base a and is de- 
noted by log a p . The most important examples of G are the multiplicative groups of a 

finite field F q , where q is either a big prime or a big power of 2 (the "classical" DLP), or 
groups emerging from the theory of elliptic curves (ECDLP). 
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Basically, there are algorithms for solving DLP that work in any group G and others, 
which take advantage of special features of the given group. 

A prominent member of the first category is the so-called baby-step giant-step algo- 
rithm by Shanks. Using the notations above, this is an algorithm that takes O(Vn) group 
operations to compute the solution x of a x =p by carrying out the following steps: 

1. Set n <- ceiling( a/Ti ) and i <-0. 

2. Compute the values a j , j=0,l,...,n-l and store them in a list. Furthermore, set 

CL ^ — CL 

3. Compare p with all elements in the list above. If p = a J for some j, then return 
x=in+j. 

4. Set p <- a p, i <- i+1 and go to step 3. 


In order to implement a simple example let's specify the group G as the prime residue 
class group mod p for some prime p and let a be an element of order n for this group. 
Now a program that computes discrete logarithms log a p in this group could look like 

this. (Note that the fourth parameter n can be omitted, if it is unknown or a is chosen 
to be a primitive root of G, as in the example below. In both cases n will be set p-1.) 

BSGS(a, p, p, n := 0, i_ := 0, j_, 1_, n_) := 

Prog 

If n = 0 
n :z p - 1 
n :z CEILING G/n) 

1_ :z ITERATES (MOD(a ■ , p) , x_, 1, n - 1) 

a := POWER_MOD(a , -n, p) 

Loop 

j_ :z POSITION^, 1_) 

If NUMBER? (j_) 

RETURN i_.n + j_ - 1 
p :z MODCccp, p) 

i_ :+ 1 


10 

NEXT_PRIME (RANDOM (10 )) = 3397745 S33 

PRIMITIVE_R00T (3397745S33) = 5 


BSGS(5, 123456789, 3397745833) = 62032454 



(11.2s) 
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Another method of this general flavour is Pollard's p -method for solving the DLP. If 
you are familiar with Pollard's p -method for the integer factoring problem, then you 
already know the general idea behind it. As for the details, have a look at the following 
program. (Note that unlike the first program the exact order n of a must be known and 
cannot be omitted here!) 

pol lard_p(o< , p, p, n, a_ := 0, b_ := 0, c_ := 0, d_ := 0, x_ := 1, y_ := 1) := 

Loop 

:= M0D(a_ + [a^, 0, l]4.(M0D(x_, 3) + 1) , n) 

b_ := MOD(b_ + [b_, 1, 0].j.(MOD(x_, 3) + 1), n) 

x_ := MOD(x_-[x_, p, tx].j. (MOD(x_, 3) + 1), p) 
c_ := MOD(c_ + [c_, 0, l!U (MOD(y_, 3) + 1), n) 

d_ := MOD(d_ + [d_, 1, 0].j. (MOD(y_, 3) +1), n) 

y_ := M0D(y_-[y_, p, ffl]|(M0D(y_, 3) + 1) , p) 
c_ := M0D(c_ + [c_, 0, l]|(M0D(y_, 3) + 1), n) 
d_ := M0D(d_ + [d_, 1, 0]|(MOD(y_, 3) +1), n) 
y_ := M0D(y_- [y_, p, a]j, (M0D(y_, 3) + 1) , p) 

If x_ = y_ 

Prog 

x_ := S0LVE_M0D((b_ - d_)-x - c_ + a_, x, n) 

Loop 

If MOD(a A FIRST(x_), p) = M0D(p, p) 

RETURN FIRST (x_) 
x := REST (x_) 


pol lard_p(5 , 123456789, 3397745833, 3397745832) = 62032454 (5.89s) 

Originally, this algorithm was only designed for the case, where n is a prime, although 
our adaption works also in the case, where n is composite. Nevertheless, in the latter 
case you should rather use the algorithm by Pohlig-Hellman, which is particularly good, 
if n is very "smooth" (see example below)! 

Pohl igHel lman(o , (i, p, n, p_, q_, x_) := 

Prog 

If PRIME?(n) 

RETURN pol lard_p(f> , (5 , p, n) 

If PRIME_POWER?(n) 

Prog 

p_ := FIRST (FIRST (FACTORS (n))) 

X- ;= pollard_p(MOD(a A (n/p_), p) , MOD(p A (n/p_) , p), p, p_) 
p := M0D(p ■ P0WER_M0D(a , -x_, p) , p) 

RETURN x_ + p_'PohligHellman(MOD(at A p_, p) , p, p, n/p_) 
q_ := VECTOR Cu_^l A u_i 2, u_, FACTORS (n)) 

CRT (VECTOR (Poh 1 i gHe 1 lman (MOD(ot A (n/ r_) , p) , M0D(p A (n/r_) , p) , p, r_) , r_, q_) , q_) 


Pohl igHel lman(5 , 123456739, 3397745333, 3397745332) = 62032454 (0.02s !!) 

3 

FACTOR (3397745 832) = 2 ■ 3 -13 -19 ■ 307 -1367 

Believe me, I would love to go into details, but there is simply no space left! If you are 
curious though, you will find them in http://www.cacr.math.uwaterloo.ca/hac/ anyway. 


